---
title: "peace_under_the_influence"
output: html_document
date: "2025-06-06"
---

```{r setup, include=FALSE}
#####libraries and data#####
library(stargazer)
library(survival)
library(tidyverse)
library(survminer)
library(MatchIt)
library(cobalt)
library(rnaturalearth)
library(rnaturalearthdata)
library(countrycode)
library(sf)
library(lmtest)
library(sandwich)
ps = read.csv("ps.csv") #### data for testing Hypothesis 1
sp = read.csv("sp.csv") #### data for testing Hypothesis 2

plot_data = read.csv("plot_data.csv")
world = st_read("world.shp")
```



```{r}
### Table 1 ###
cox1 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_both + intensity + duration + rebstrength_num + terrcont_num + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)


cox2 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_both + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)


cox3 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_both + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)

stargazer(cox1, cox2, cox3, apply.coef = exp,  p.auto = FALSE, 
          t.auto = FALSE, covariate.labels = c(
    "Drug",
    "Conflict Intensity",
    "Conflict Duration",
    "Rebel Strength X Time",
    "Territorial Control",
    "GDP per Capita",
    "Democratizing Move"
  ))


#### Table A3 ####
cox4 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_extortion + intensity + duration + rebstrength_num + terrcont_num + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)


cox5 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_extortion + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)


cox6 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_extortion + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)

stargazer(cox4, cox5, cox6, apply.coef = exp,  p.auto = FALSE, 
          t.auto = FALSE, covariate.labels = c(
    "Drug Extortion",
    "Conflict Intensity",
    "Conflict Duration",
    "Rebel Strength X Time",
    "Territorial Control",
    "GDP per Capita",
    "Democratizing Move"
  ))

#### Table A4 ####
cox7 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_smuggling + intensity + duration + tt(rebstrength_num) + terrcont_num + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)


cox8 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_smuggling + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)


cox9 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_smuggling + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)

stargazer(cox7, cox8, cox9, apply.coef = exp,  p.auto = FALSE, 
          t.auto = FALSE, covariate.labels = c(
    "Drug Smuggling",
    "Conflict Intensity",
    "Conflict Duration",
    "Rebel Strength X Time",
    "Territorial Control",
    "GDP per Capita",
    "Democratizing Move"
  ))


#### Table 5 ####
cox_oil = coxph(Surv(start, stop, first_peace_process_start) ~ 
       oil + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)

cox_gem = coxph(Surv(start, stop, first_peace_process_start) ~ 
       gem + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)


cox_gold = coxph(Surv(start, stop, first_peace_process_start) ~ 
      gold + intensity + duration + tt(rebstrength_num) + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = ps,  tt = function(x, t, ...) x * t)


stargazer(cox_oil, cox_gem, cox_gold, apply.coef = exp,  p.auto = FALSE, 
          t.auto = FALSE, covariate.labels = c(
    "Oil",
    "Gem",
    "Gold",
    "Conflict Intensity",
    "Conflict Duration",
    "Rebel Strength X Time",
    "Territorial Control",
    "GDP per Capita",
    "Democratizing Move"
  ))

```


```{r}
logit1 = glm(peace_process_start ~ drug_both + intensity + duration + rebstrength_num + terrcont_num, data = ps, family = "binomial")
  
logit2 = glm(peace_process_start ~ drug_both + intensity + duration + rebstrength_num + terrcont_num + gdppc, data = ps, family = "binomial")


logit3 = glm(peace_process_start ~ drug_both + intensity + duration + rebstrength_num + terrcont_num + gdppc + demo_move, data = ps, family = "binomial")



logit1_robust = sqrt(diag(vcovCL(logit1,  cluster = ps$dyad_id, type = "HC1")))

logit2_robust = sqrt(diag(vcovCL(logit2,  cluster = ps$dyad_id, type = "HC1")))
logit3_robust = sqrt(diag(vcovCL(logit3,  cluster = ps$dyad_id, type = "HC1")))


stargazer(logit1, logit2, logit3, 
  type = "latex",
  se = list(logit1_robust, logit2_robust, logit3_robust),
  covariate.labels = c(
   "Drug",
    "Conflict Intensity",
    "Conflict Duration",
    "Rebel Strength",
    "Territorial Control",
    "GDP per Capita",
    "Democratizing Move"
  )
)
```

```{r}
#### Figure 3 ####
km_obj = Surv(time = ps$stop, event = ps$event)

ggsurvplot(
  survfit(km_obj ~ drug_both, data = ps),
  data = ps,
  conf.int = TRUE,
  pval = TRUE,
  xlab = "Time (Years)",
  ylab = "Survival Probability",
  legend.title = "",
  legend.labs = c("No Drug", "Drug"),
  ggtheme = theme_minimal(),
  ylim = c(0.5, 1),
  pval.coord = c(1, 0.7),
  xlim = c(0, 10)
)
```


```{r}
#### Table 2 ####
splinter1 = glm(splinter_emerge ~ drug_both + rebel_age + terrcont_num +  rebstrength_num +ties_dummy, family = "binomial", data = sp)
splinter2 = glm(splinter_emerge ~ drug_both + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity  + multiparty_conflict, family = "binomial", data = sp)

splinter3 = glm(splinter_emerge ~ drug_both + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial, family = "binomial", data = sp)

splinter_robust1 = sqrt(diag(vcovCL(splinter1, cluster = sp$dyad_id, type = "HC1")))
splinter_robust2 = sqrt(diag(vcovCL(splinter2, cluster = sp$dyad_id, type = "HC1")))
splinter_robust3 = sqrt(diag(vcovCL(splinter3, cluster = sp$dyad_id, type = "HC1")))

stargazer(splinter2, se = list(splinter_robust2))
stargazer(splinter1, splinter2, splinter3,
  type = "latex",  # or "text"
  se = list(splinter_robust1, splinter_robust2, splinter_robust3),
  covariate.labels = c(
    "Drug",
    "Rebel Age",
    "Territorial Control",
    "Rebel Strength",
    "Social Ties",
    "Conflict Intensity",
    "Multiparty Conflict",
    "Territorial Conflict"
  ),
  column.labels = c("Model 1", "Model 2", "Model 3")
)


#### Table A5 ####
se1 = glm(splinter_emerge ~ drug_extortion + rebel_age + terrcont_num +  rebstrength_num +ties_dummy, family = "binomial", data = sp)
se2 = glm(splinter_emerge ~ drug_extortion + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity  + multiparty_conflict, family = "binomial", data = sp)

se3 = glm(splinter_emerge ~ drug_extortion + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial, family = "binomial", data = sp)

se_robust1 = sqrt(diag(vcovCL(se1, cluster = sp$dyad_id, type = "HC1")))
se_robust2 = sqrt(diag(vcovCL(se2, cluster = sp$dyad_id, type = "HC1")))
se_robust3 = sqrt(diag(vcovCL(se3, cluster = sp$dyad_id, type = "HC1")))


stargazer(se1, se2, se3,
  type = "latex",  
  se = list(se_robust1, se_robust2, se_robust3),
  covariate.labels = c(
    "Drug Extortion",
    "Rebel Age",
    "Territorial Control",
    "Rebel Strength",
    "Social Ties",
    "Conflict Intensity",
    "Multiparty Conflict",
    "Territorial Conflict"
  ),
  column.labels = c("Model 1", "Model 2", "Model 3")
)


#### Table A6 ####
sm1 = glm(splinter_emerge ~ drug_smuggling + rebel_age + terrcont_num +  rebstrength_num +ties_dummy, family = "binomial", data = sp)
sm2 = glm(splinter_emerge ~ drug_smuggling + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity  + multiparty_conflict, family = "binomial", data = sp)

sm3 = glm(splinter_emerge ~ drug_smuggling + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial, family = "binomial", data = sp)
sm_robust1 = sqrt(diag(vcovCL(sm1, cluster = sp$dyad_id, type = "HC1")))
sm_robust2 = sqrt(diag(vcovCL(sm2, cluster = sp$dyad_id, type = "HC1")))
sm_robust3 = sqrt(diag(vcovCL(sm3, cluster = sp$dyad_id, type = "HC1")))


stargazer(sm1, sm2, sm3,
  type = "latex",  
  se = list(sm_robust1, sm_robust2, sm_robust3),
  covariate.labels = c(
    "Drug Smuggling",
    "Rebel Age",
    "Territorial Control",
    "Rebel Strength",
    "Social Ties",
    "Conflict Intensity",
    "Multiparty Conflict",
    "Territorial Conflict"
  ),
  column.labels = c("Model 1", "Model 2", "Model 3")
)


#### Table 6 ####
oil_splinter = glm(splinter_emerge ~ oil + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial, family = "binomial", data = sp)
gem_splinter = glm(splinter_emerge ~ gem + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial, family = "binomial", data = sp)
gold_splinter = glm(splinter_emerge ~ gold + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial, family = "binomial", data = sp)


oil_robust  = sqrt(diag(vcovCL(oil_splinter,  cluster = sp$dyad_id, type = "HC1")))
gem_robust  = sqrt(diag(vcovCL(gem_splinter,  cluster = sp$dyad_id, type = "HC1")))
gold_robust = sqrt(diag(vcovCL(gold_splinter, cluster = sp$dyad_id, type = "HC1")))

stargazer(oil_splinter, gem_splinter, gold_splinter,
  type = "latex",  
  se = list(oil_robust, gem_robust, gold_robust),
  covariate.labels = c(
    "Oil",
    "Gem",
    "Gold",
    "Rebel Age",
    "Territorial Control",
    "Rebel Strength",
    "Social Ties",
    "Conflict Intensity",
    "Multiparty Conflict",
    "Territorial Conflict"
  ),
  column.labels = c("Oil", "Gem", "Gold")
)

```


```{r}
#### Figure 4 ####
get_mode = function(x) {
  ux = na.omit(unique(x))
  ux[which.max(tabulate(match(x, ux)))]
}

newsplinter = sp %>%
  summarise(
    rebel_age = get_mode(rebel_age),
    terrcont_num = median(terrcont_num, na.rm = TRUE),
    rebstrength_num = median(rebstrength_num, na.rm = TRUE),
    ties_dummy = median(ties_dummy, na.rm = TRUE),
    intensity = median(intensity, na.rm = TRUE),
    multiparty_conflict = median(multiparty_conflict, na.rm = TRUE),
    territorial = mean(territorial, na.rm = TRUE)
  ) %>%
  slice(rep(1, 2)) %>%
  mutate(drug_both = c(0, 1))


newsplinter$pred = predict(splinter3, newdata = newsplinter, type = "response") * 100


library(ggplot2)

ggplot(newsplinter, aes(x = factor(drug_both, labels = c("No Drug", "Drug")), y = pred)) +
  geom_col(width = 0.5, fill = "firebrick") +
  geom_text(aes(label = sprintf("%.2f%%", pred)), vjust = -0.5, size = 4) +
  labs(
    x = "Drug Activities",
    y = "Predicted Probability of Splintering (%)",
    title = ""
  ) +
  ylim(0, max(newsplinter$pred) + 1) +
  theme_minimal()


```


```{r}
##### Matching #####

#### Table 3 ####
m_ps3 = matchit(drug_both ~ intensity + duration + rebstrength_num + terrcont_num + gdppc + demo_move,
                 data = ps, method = "full", mahvars = c("intensity", "duration", "gdppc"), caliper = 0.9)

m_ps4 = matchit(drug_both ~ intensity + duration + rebstrength_num + terrcont_num + gdppc + demo_move,
                 data = ps, method = "full", mahvars = c("terrcont_num", "duration", "gdppc"), caliper = 0.55)

matched_ps3 = match.data(m_ps3)
matched_ps4 = match.data(m_ps4)

cox_matching1 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_both + intensity + duration + rebstrength_num + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = matched_ps3,  tt = function(x, t, ...) x * t, weights = weights)



cox_m1 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_both + tt(intensity) + duration + rebstrength_num + tt(terrcont_num) + gdppc + demo_move + cluster(dyad_id), data = matched_ps3,  tt = function(x, t, ...) x * t, weights = weights)



cox_matching2 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_both + intensity + duration + rebstrength_num + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = matched_ps4,  tt = function(x, t, ...) x * t, weights = weights)


cox_m2 = coxph(Surv(start, stop, first_peace_process_start) ~ 
       drug_both + intensity + tt(duration) + rebstrength_num + terrcont_num + gdppc + demo_move + cluster(dyad_id), data = matched_ps4,  tt = function(x, t, ...) x * t, weights = weights)

stargazer(cox_m1, cox_m2, apply.coef = exp,  p.auto = FALSE, 
          t.auto = FALSE, covariate.labels = c(
    "Drug",
    "Conflict Intensity X Time",
    "Conflict Duration",
    "Conflict Intensity",
    "Conflict Duration X Time",
    "Rebel Strength",
    "Territorial Control X Time",
    "Territorial Control",
    "GDP per Capita",
    "Democratizing Move"
  ))
```


```{r}
#### Figure A1
love.plot(m_ps3, 
          stats = "mean.diffs", 
          threshold = 0.1, 
          var.order = "unadjusted",
          abs = TRUE, 
          line = TRUE,
          colors = c("black", "firebrick"),
          shapes = c("circle", "triangle"), var.names = c(
            "distance" = "Distance",
    "intensity" = "Conflict Intensity",
            "rebstrength_num" = "Rebel Strength",
            "terrcont_num" = "Territorial Control",
            "demo_move" = "Democratizing Move",
            "duration" = "Conflict Duration",
            "gdppc" = "GDP per Capita"
          ))

#### Figure A2
love.plot(m_ps4, 
          stats = "mean.diffs", 
          threshold = 0.1, 
          var.order = "unadjusted",
          abs = TRUE, 
          line = TRUE,
          colors = c("black", "firebrick"),
          shapes = c("circle", "triangle"), var.names = c(
            "distance" = "Distance",
    "intensity" = "Conflict Intensity",
            "rebstrength_num" = "Rebel Strength",
            "terrcont_num" = "Territorial Control",
            "demo_move" = "Democratizing Move",
            "duration" = "Conflict Duration",
            "gdppc" = "GDP per Capita"
          ))



```

```{r}
#### Table 4
m_splinter3 = matchit(drug_both ~ rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial,
                 data = sp, method = "full", caliper = 0.1, mahvars = c("multiparty_conflict", "rebel_age", "terrcont_num"))



m_splinter4 = matchit(drug_both ~ rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial,
                 data = sp, method = "full", caliper = 0.5, mahvars = c("rebel_age"))


matched_splinter3 = match.data(m_splinter3)
matched_splinter4 = match.data(m_splinter4)

splinter3_matched = glm(splinter_emerge ~ drug_both + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial, weights = weights,
                   family = "binomial", data = matched_splinter3)

splinter4_matched = glm(splinter_emerge ~ drug_both + rebel_age + terrcont_num +  rebstrength_num + ties_dummy + intensity   + multiparty_conflict + territorial, weights = weights,
                   family = "binomial", data = matched_splinter4)

splinter3_robust <- sqrt(diag(vcovCL(splinter3_matched, cluster = matched_splinter3$dyad_id, type = "HC1")))
splinter4_robust <- sqrt(diag(vcovCL(splinter4_matched, cluster = matched_splinter4$dyad_id, type = "HC1")))


stargazer(splinter3_matched, splinter4_matched,
  type = "latex", 
  se = list(splinter3_robust, splinter4_robust),
  covariate.labels = c(
    "Drug",
    "Rebel Age",
    "Territorial Control",
    "Rebel Strength",
    "Social Ties",
    "Conflict Intensity",
    "Multiparty Conflict",
    "Territorial Conflict"
  ),
  column.labels = c("Model 1", "Model 2")
)
```


```{r}
#### Figure A3
love.plot(m_splinter3, 
          stats = "mean.diffs", 
          threshold = 0.1, 
          var.order = "unadjusted",
          abs = TRUE, 
          line = TRUE,
          colors = c("black", "firebrick"),
          shapes = c("circle", "triangle"), var.names = c(
            "intensity" = "Conflict Intensity",
            "territorial" = "Territorial Conflict",
            "multiparty_conflict" = "Multiparty Conflict",
            "rebstrength_num" = "Rebel Strength",
            "ties_dummy" = "Social Ties",
            "terrcont_num" = "Territorial Control",
            "demo_move" = "Democratizing Move",
            "distance" = "Distance",
            "rebel_age" = "Rebel Age"
          ))

#### Figure A4
love.plot(m_splinter4, 
          stats = "mean.diffs", 
          threshold = 0.1, 
          var.order = "unadjusted",
          abs = TRUE, 
          line = TRUE,
          colors = c("black", "firebrick"),
          shapes = c("circle", "triangle"), var.names = c(
            "intensity" = "Conflict Intensity",
            "territorial" = "Territorial Conflict",
            "multiparty_conflict" = "Multiparty Conflict",
            "rebstrength_num" = "Rebel Strength",
            "ties_dummy" = "Social Ties",
            "terrcont_num" = "Territorial Control",
            "demo_move" = "Democratizing Move",
            "distance" = "Distance",
            "rebel_age" = "Rebel Age"
          ))
```

```{r}
#### Figure 1

ggplot(data = world) +
  geom_sf(aes(fill = drug_group), color = "gray80") +
  scale_fill_manual(values = c("Yes" = "firebrick", "No" = "gray90"), name = "Rebel Groups Financing via Illicit Drugs") +
  theme_minimal() +
  labs(title = "",
       subtitle = "",
       caption = "") +
  theme(legend.position = "bottom")

#### Figure 2
plot_data$drug_use = factor(plot_data$drug_use, levels = c("No Drug", "Drug"))
ggplot(plot_data, aes(x = factor(country), y = n_groups, fill = drug_use)) +
  geom_col(position = "stack") +
  scale_fill_manual(values = c("No Drug" = "gray70", "Drug" = "firebrick")) +
  labs(x = "", y = "Number of Rebel Groups",
       title = "",
       fill = "Financing via Illicit Drugs") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 5))
```




